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ABSTRACT 

It is known that the stability of finite difference models of hyperbolic 
initial boundary value problems is connected with the propagation and 
reflection of parasitic waves. Here the waves point of view is applied to 
models containing two boundaries or interfaces, \diere repeated reflection of 
trapped wave packets is a potential new source of instability. Our analysis 
accounts for various known instability phenomena in a unified way and leads to 
several new results, three of which are as follows. (1) Dissipativity does 
not ensure stability when three or more formulas are concatenated at a 
boundary or internal interface. (2) Algebraic "GKS instabilities" can be 
converted by a second boundary to exponential instabilities only when an 
infinite numerical reflection coefficient is present. (3) "GKS-stability" and 
"P“Stability" can be established in certain problems by showing that all 
numerical reflection coefficients have modulus less than 1. 
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Ufotation 


j* n 
h, k 
^ = k/h 
u(x,t), 

N 

C(C,w) 

z 

K, [X 
A, r 

R = NJ?., L = Nr 

Q, Q 
A(z) 
const 


space, time index 
space, time step size 
mesh ratio 

continuous, discrete solution vector 
dimension of u, v 
wave number, frequency 
group velocity 

temporal amplification factor 

rightgoing, leftgoing spatial amplification factor 
no. of pts. to the left, right of center in stencil 
no. of rightgoing, leftgoing numerical wave modes 
difference model for i.v.p., i.b.v.p. 
reflection matrix function at boundary 
positive constant, different each time 
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0. INTRODUCTION 


Hyperbolic systems of partial differential equations admit solutions which 
behave locally like waves moving along characteristics. When such a system is 
modeled numerically by finite differences or finite elements, the result is a 
dispersive medium that may admit additional parasitic wave modes with 
wavelengths on the scale of the discretization. Energy associated with these 
parasitic waves travels at a group velocity that is unrelated to the 
characteristics of the original system [25], [30]. However the behavior of 

such waves has a decisive effect on stability. 

For finite difference models of linear hyperbolic problems with a single 
spatial boundary, a stability theory was developed around 1970 by Kreiss , 
Osher, and others [10], [18]. In earlier papers we have shown that this 

theory can be naturally stated in terms of dispersive wave propagation [26], 
[27]. To summarize: if a boundary with homogeneous boundary conditions can 

emit a radiated wave in the absence of any incident waves, i.e., a wave with 
group velocity pointing into the interior of the domain, then it is 
unstable. If it has an infinite reflection coefficient for waves at some 
frequency, a stronger condition, then it is more severely unstable. 

This paper applies wave propagation ideas to investigate stability for 
one-dimensional linear finite difference models with two or more boundaries or 
internal Interfaces. The most basic example of such a model is a discrete 
approximation to an equation whose spatial domain is an interval such as 
[0,1]. Another example is a model of a problem featuring discontinuous 
coefficients, e.g. wave propagation in a discontinuously stratified medium 
[4]. A third is a numerical scheme employing local mesh refinement to improve 
in which various interfaces between fine and coarse grids will be 


accuracy , 


present [2] , [5] • A fourth is any model with a composite numerical boundary 

or interface, such as a fourth-order difference model on [0,®) that has a 
five-point stencil, and which therefore requires one numerical boundary 
condition at j = 0 and another at j = 1 [15]. Such a model can be viewed 
as containing two interfaces separated by a single grid point, and we will 
show that this view may be useful for stability analysis. 

Any multi— interface model potentially admits trapped numerical waves that 
reflect back and forth repeatedly from one interface to another. If the 
reflections cause amplification, this can lead to unbounded growth of 
numerical solutions. The factors that control this are: magnitude of the 

reflection coefficients, which is the source of growth; dissipation of waves 
as they travel between interfaces, which is a source of attenuation; and 
travel time between interfaces, which determines how frequently any reflection 
circuit that causes growth is repeated. All of the arguments of this paper 
consist of working out consequences of various combinations of these factors 
that may be of practical interest. 

In particular we investigate two kinds of stability. First, "stability" 
or "Lax-Richtmyer stability" refers to the usual Lax-Richtmyer definition for 
time-dependent finite difference models, or to variants thereof such as "GKS- 
stability" (Defn. 3.3 in the well-known paper by Gustafsson, Krelss, and 
Sundstrom [10]). A difference model that is stable in this sense may admit 
solutions that grow with time, provided that the growth does not get worse as 
the mesh is refined. This is what is needed to ensure convergence as the mesh 
size approaches zero to the correct solution of the time-dependent 
differential equation, for each fixed time t. On the other hand to be "P- 
stable" [1], a model must admit no growing solutions at all. (See Section 3 
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for the precise definition. Such a model is also sometimes called '’time- 

stable** [29].) Although the theory here is not as well developed, such a 
condition is needed if a time-dependent difference model is to be used to 
obtain steady-state solutions, as is common in computational aerodynamics. As 
a rule of thumb, we will show that P-instability is very often associated with 
reflection coefficients greater than 1 in magnitude, and Lax-Richtmyer 
instability with reflections coefficients that are infinite. 

Section 1 reviews stability theory for one-boundary difference models 
(Prop. 1). Section 2 investigates interfaces separated by a fixed number of 
grid points Aj as the mesh is refined, as in the fourth-order boundary 
condition mentioned above. Here the travel times go to zero with the mesh 
spacing, with the effect that finite reflection coefficients greater than 1 
can cause catastrophic unstable growth, regardless of dissipation (Props. 3, 
4, 4^). Conversely, reflection coefficients smaller than 1 ensure stability 
(Prop. 5). Section 3 considers interfaces separated by a fixed distance Ax 
as the mesh is refined, as in the problem on [0,1] mentioned above. Here the 
travel times are independent of mesh spacing, so large finite reflection 
coefficients can cause P-instability (Prop. 7) but not Lax-Richtmyer 
instability (Prop. 6). In this context multiple reflections may convert the 
weak instability of a single interface to a catastrophic instability (Prop. 
8), but only if the unstable interface is of the sort with an infinite 
reflection coefficient (Prop. 9). In contrast to the fixed-Aj case, 
dissipation now prevents all kinds of solution growth (Prop. 10), and so again 
do reflection coefficients less than 1 (Prop. 11). 

For convenience of reference, here is a list of the explicit examples 
presented here to illustrate various points. The symbol A indicates the 
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modulus of a reflection coefficient, and the solution norm at time step 

n. These quantities will be made more precise later on. 

1* Algebraically unstable one~boundary model (one boundary, A = , 

S const^l. 
n ' 

2. Exponentially unstable concatenation of three stable dissipative formulas 
(fixed-Aj, A > 1, S ~ constn). 

3. P-stability guaranteed by reflection coefficients less than 1 

(fixed-Aj or Ax, A < 1, ~ const). 

P-instability caused by reflection coefficients greater than 1 
(fixed-Ax, A > 1, S ~ const*”). 

5. Exponential instability caused by interaction of two algebraically 
unstable boundaries [fixed-Ax, A = «>, ~ (Aj)*^®”^*" *”) . 

The reader may be disappointed at the artificiality of some of these 
examples, especially (2) and (3), and he may wonder how helpful wave 
reflection ideas can be in practice for the design of difference schemes. A 
full answer to this question will have to await further experience. 
Nevertheless there is no doubt that the instability mechanisms described here 
are real and deserve to be understood. At present, virtually no difference 
models have been shown to be stable of the sort that contain multiple 
interfaces at various points, such as might appear in adaptive mesh 
refinement. Perhaps the ideas here, such as Prop. 5, can help bring about a 
change in this situation. 

Much of the material in this paper can be found in Section 6 of the 
author s Ph.D. dissertation [24]. For some numerical illustrations, see [28]. 
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For a different analysis of stability for two-boundary problems that is 
closely related to the present one, see the report [8] by Giles and Thompkins , 
which is mainly concerned with P-stability« Giles and Thompkins consider 
parasitic wave propagation for models with variable as well as constant 
coefficients . 

The phenomenon of instability caused by trapped wave packets can also 
occur in two-dimensional problems when the domain contains a corner. Osher 
has given examples of hyperbolic systems (not difference models) in corners 
that are ill-posed because of trapped waves [19] , while Sarason and Smoller 
have shown that for a 2^2 strictly hyperbolic system such as the second-order 
wave equation, this cannot happen [21] . But trapped numerical waves may 
render a finite difference model of even the latter sort unstable. The 
principles involved are precisely those of this paper, but we will discuss 
corners elsewhere. 

The reader interested in getting to the main ideas quickly may find it 
possible to turn directly to Section 2. 

1. REVIEW OF WAVE PROPAGATION AND STABILITY FOR ONE-BOUNDARY DIFFERENCE IfflDELS 

Consider a linear first-order hyperbolic system of differential equations 

u = Au (!• 1) 


with initial data 


u(x,0) = f (x) • 


( 1 . 2 ). 


Here u(x,t) and f(x) are N-vectors, A is a constant NXN matrix, and the 
spatial domain is ]R. The statement that (1»1) is hyperbolic means that A 
has real eigenvalues Kv<N, and a complete set of associated 

eigenvectors • It follows that if C £ is an arbitrary wave number , 

then (1.1) admits N linearly independent solutions of the form 
u(x,t) = Uexp(i(cot + 5x) ) , namely the waves 

i(w^(?)tHx) 

u(x,t) = e ( 1 . 3 ) 

with C. 03 ^ is called the frequency of (1.3), and the N-valued 

linear function w = is the dispersion relation for (1.1). Each wave 

(1.3) propagates uniformly with no change in shape at the velocity -p^, hence 
leftward or rightward depending on whether p^ is positive or negative, 
respectively . 

Since the vectors span any f e L^( :b?^) can be written as a 

Fourier Integral with respect to ^ of waves (1.3). It follows by Parseval's 
theorem that Hu(*,t)ll is constant with respect to t; a fortiori, for any 
fixed t one has 

flu(*,t)H < const llfll, (1.4) 

which is to say, (1.1) - (1.2) is well posed in L^. Related well-posedness 
bounds continue to hold under reasonable assumptions if (1.1) is given a 
zeroth-order term Bu, an inhomogeneous term F(x,t), or variable 

coefficients, although in these circumstances some growth at a bounded rate 
in t must be permitted. For simplicity we will ignore these possibilities. 
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Let u be approximated by a vector grid function v? = v(jh,nk) 

» u(jh,nk), where k is a time step and h is a space step, {v”} will be 
determined iteratively as the solution of an s+2— level finite difference 
formula 



n+1 

V 



n -0 

V 


(1.5) 


where each Q is a spatial finite difference operator with constant matrix 
coefficients of dimension NxN* Let Q be a name for (1*5). As with the 
differential equation, one can show that Q admits solutions 


^ ^^i((ot+5x)^ X = jh, t = nk, V e (1.6) 

For each 5 e IR, in fact, it permits in general not N but (s+l)N distinct 
values of a>, whose relation to C constitutes the dispersion relation for 
(1.5). These values depend nonlinearly on 5, and they are not necessarily 
real. A solution with 5 ^ ]R and Imu) > 0 decays with t, but a solution 
with 5 G IR and Inw < 0 grows at the rate = const”, and if Q 

admits a solution of this kind, it is unstable. On the other hand if there 
are no such growing modes, and if any modes with ^ ,to e are nondefective 

in a sense we will not go into, then Q is stable. Thus stability for a 
constant coefficient finite difference model on an unbounded domain can be 
investigated by a fairly straightforward process of Fourier analysis. For 
details, see [20]. 

Let Q be stable and admit a solution (1.6) with 5 e IR. It can be 
shown that the dispersion relation for (1.5) determines a function to = to(l') 
for 5 in a neighborhood of [27, Lemma 3.2], and that the energy 
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associated with the wave (1.6) propagates at the group velocity 


C(C,w) 


d(jJ 

dC 


(C,W) . 


(1.7) 


The precise meaning of this statement is asymptotic: if Q is given Initial 

data 


fj = <Kx-Ct)Ve^^“*^‘*'^^^ 


0 < a < s, t = ok, X = jh 


for some smooth function (|), then the solution at a later time will be 


v^ » <Kx-Ct)Ve^^‘^^+^^> 


with the equality becoming more exact as (t> is made smoother. See for 
example Lemma 5.1 of [27] . 

Example 1. As an example, consider the leap frog (LF) model 







const 


( 1 . 8 ) 


of u^ = u^. By inserting v(x,t) = , one finds that the dispersion 

relation is 


sin wk = X sin 5h. 


(1.9) 


Differentiation leads to the group velocity formula 


ca,oi) = - 


cos Ch 
cos tok 


( 1 . 10 ) 
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Thus a well-resolved wave, i.e. one with Ch, lok “ 0, propagates under LF with 
group velocity C “ -1* On the other hand LF also admits many waves with Ch 
or o)k not small. The extreme cases are the "parasitic" solutions 
(C,w) = (Tt/h.O), (0,i:/k), and (n/h,ii/k) , which by (1.10) have group veloci- 
ties +1, +1, and -1, respectively. For the first two of these, the sign of 
C reveals that energy propagates in the physically wrong direction. In fact 
for each sufficiently small frequency u e E, (1.9) gives two distinct wave 
numbers ^ in the fundamental range [-ii/h,x/h] , and by (1.10), one of the 
corresponding waves propagates leftwards and the other propagates rightwards. 

See [25] or [30] for illustrations. // 

Returning to the general model Q of (1.5), let us change the notation 

and rewrite (1.6) in the more convenient form 

z” K,z e (t-{0} (1.11) 

where < = e^^'^ and z = e^“^. (For full generality one must permit a 

further multiplicative factor j to represent certain defective modes. Such 

modes are rarely important in practice, however, so in all of what follows we 
assume 6=0, although the results remain valid in the general situation. 
The reader is referred to [27] for more complete details.) A solution (1.11) 
with |k1=|z| = 1 and C<0 (resp. > 0) will be called leftgoing (resp . 
rightgoing) . For obvious geometric reasons it also makes sense to say that a 
solution with |z| > 1 Is leftgoing If |k| > 1 and rightgoing if |k| < 1. 
It can be shown under reasonable assumptions (see [10]) that for any z with 
|z| > 1, Q admits a family of R = Ni linearly independent rightgoing and 
L = Nr linearly independent leftgoing solutions (1.11), where A and r 
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denote the numbers of grid points to the left and right of center, 
respectively, covered by the stencil of Q [27, Sec. 3]. Therefore, the 
general solution to (1.5) of the form <1)^ is a linear combination 

n „ R+L 

V = z I a V . 

3 , m m m 

m=l 

If we relabel a, V, and k by P, W, and \x for leftgoing components, this 
becomes 


n n V i n V 

V = z 2, “ V K-' + z y P W 

J m m m m : 


m m m 


( 1 . 12 ) 


(RIGHTGOING) 


(LEFTGOING) 


We emphasize that the leftgoing and rightgoing waves in this sum have very 
little to do with the waves admitted by the original equation (1.1). 

Let a lefthand boundary be introduced at x = 0, so that the spatial 
domain becomes :r+ and j is restricted to j > 0. Now (1.1) must be 

supplemented by as many additional scalar boundary conditions as there are 
inflowing characteristics at x = 0, and if this done in the natural way, 
well-posedness is guaranteed [13]. But we pass over these details and 
consider the finite difference model. Regardless of the characteristics of 
(1.1), (1.5) will have to be supplemented by R = NJl additional boundary 

conditions, one for each rightgoing numerical solution component. For 
simplicity we take these to be homogeneous and of the form 


n+1 


1 “2 

I I 

cr=-l i=0 


n-o 


0 < j < A-1 (1.13) 
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for some integers and M 2 and N^N matrices Q ^ name for 

the initial boundary value problem model (1.5), (1.13). 

In practice it can easily happen that Q is unstable. A theory of such 
instability was developed a decade ago by Kreiss, Osher, and others, and 

described at length in the well-known paper [10] by Gustafsson, Kreiss, and 
Sundstrom — henceforth "GKS". In [26] and [27] the Kreiss/Osher theory has 
been given the following interpretation. If |z| > 1 is fixed, then the 

general superposition (1.12) of leftgoing and rightgoing waves does not 
satisfy (1.13), and hence is not a solution to Q. Instead, (1.13) can be 

thought of as a set of R = Nl reflection conditions relating rightgoing to 

leftgoing waves at the boundary. These conditions are obtained by substitu- 
ting (1.12) in (1.13) and then collecting terms in and so that one 

gets 



“1 


“Pl 


• 

= D(Z) 

• 

E(z) 

• 

• 

• 

• 


_“r_ 




for some R^R matrix E(z) and R^L matrix D(z) . For most z, E(z) will be 
nonsingular, and (1.14) determines the reflected wave coefficients as a 
bounded function of the incident ones. If we write A = E ^ D, so that 

A(z) is the RXL reflection matrix for the given boundary conditions, then 
this function has the form 

a = A(z)P = [E(z)]”^ D(z)P. (1.15) 
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(This A(z) has nothing to do with the coefficient matrix of (1.1).) 
However, it may happen that for some |zq| > 1, E(zq) Is singular, and in 
this case (1.14) permits a solution consisting of rightgoing waves in the 
absence of leftgoing waves. This will cause instability. If in this 

situation A(z) is unbounded as z ^ Zq, then an infinite reflection 

^o®f_f_ic^en^ exists at and the instability will be particularly severe 


Thus the Kreiss/Osher theory leads to the following "GKS stability 
theorem" ; 


PROPOSITION 1 [10, 27] . Ihe initial boundary value problem model Q is 

GKS-uns table if and only if E(z) is singular for some |z| > 1. 
Equivalently, it is GKS-unstable if and only if for some |z| > 1 it admits a 
nonzero — solution v” = z” 4>j (1.12) consisting entirely of rightgoing wave 

components . 

Proof . See [27] for a precise statement and proof. 


The notion of GKS— stability" employed in this result is a fairly complicated 
one given as Defn. 3*3 in [10] . See [27] for a discussion of its meaning. For 
the remainder of this paper "stable" means "GKS-stable", except where 
otherwise stated. 

Example 1, co ntinued . To return to the previous example, suppose LF 
is applied on x > 0 with the numerical boundary condition 


^n+1 ^ n+1 

0 1 


(1.16) 
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In K,z notation, the dispersion relation (1.9) and group velocity (1.10) for 


LF are 


z 



_ < + l/< 

” z + 1/z ’ 


(1.17) 


and (1.16) imposes the additional condition < = 1. One sees immediately 

that the wave (k»z) = (1»-1)> i.e., v^ = (-1)^> satisfies both the interior 

formula and the boundary condition and has C > 0. Therefore by Prop. 1 the 

, . , n+1 n . 

model (1.8), (1.16) is unstable. By contrast the condition Vq = v^ is 

satisfied by no rightgoing solutions to LF, so with this boundary condition 
LF would be stable. 

This example is one of those with an infinite reflection coefficient. To 
see this, note that for each |z| > 1, (1.17) gives two values of k related 
by <2 = -I/Kj. Let these be denoted by < and p, where k is the "right- 
going" value with Re< Rez ^ 0 and |k| ^ 1, for which C ^ 0 if |k| — 1, 

and p is the "leftgoing" one with ReK Rez > 0 and |k| > 1. Then for 

this problem the superposition (1.12) takes the form v^ = ocp^ + . To 

calculate the reflection coefficient we substitute this in (1.16) and obtain 
a + p = ap + Pk, that is, p = Aa with 


This quotient becomes infinite when z=-l, k= 1, p=-l. 

The unstable behavior of this difference model is illustrated in Figs. 4.1 
-4.2 and Figs. 5.1 - 5.4 of [24] and in Figs. 3, 4 of [27]. // 
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2. TWO INTERFACES SEPARATED BY A FIXED NUMBER OF GRID POINTS Aj 

The stability result of Prop. 1 is illustrated in Fig. 1. If a set of 
numerical waves reflects at a boundary with a gain in amplitude, as in Fig. 
la, this does not constitute instability. It may force the constant in a 
discrete estimate analogous to (1.4) to be large, but it does not preclude the 
existence of such an estimate. On the other hand if the boundary can produce 
radiated energy in the presence of no incident energy at all, as in Fig. lb, 
then it is unstable. 






(a) stable 


(b) unstable 


Figure 1. 


Stable and unstable solutions 


*3 


at a lefthand boundary. 


Suppose now that Q is a model containing not a boundary but an internal 
interface of some kind separating two difference schemes Q_ and Q_|_ 

(possibly identical) • The Interface might be a complicated structure 

extending over several grid points, or it might be simply an abrupt change of 
coefficient, of difference formula, or of mesh size. It is plausible that the 
picture should change to that of Fig. 2: Q is unstable if and only if it 

permits a solution z (|) that is outgoing from the interface on both 
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sides. This conclusion can be derived rigorously from Prop. 1 by folding the 
interface problem into an initial boundary value problem for a system of 
equations of twice the original size [3!, [5] » [6], [16], [24]. 



<— 




Figure 2. 


(a) stable 


Stable and unstable solutions 


(b) 


n A 
Z <P, 


unstable 


at an Internal interface. 


Reflection equations for an internal interface analogous to (1.14) - 

(1.15) for a boundary can be obtained by the same folding idea. For each 
|zl > 1, there are R~ + L"^ linearly independent waves that may be incident 
at the interface from both sides, and L + R"*" that may be radiated. The 
full reflection equation is the linear system describing how the coefficients 
of these wave components are related. 




“l 

• 


• 

• 


• 

• 


• 

r 


a“ 

L 

= D(z) 

r" 

+ 


Pt 

a 

1 


• 

• 


• 

• 

• 


• 




R+ _ 
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where E and D are matrix functions of dimensions (L“ + R'*’) x (l“ + r+) 

and (L + R ) x (r + L ) (cf. (1.14)). However, in this paper we will only 
need the response of an interface on one side to a wave incident on that side. 
The corresponding reflection equation is the projection of (2.1) onto a one- 
sided domain and range. In the case of incidence on the right, for example, 
it has the form 



where E is R+ x R+ and D is R+ x l+. when E(z) is nonsingular, 

(2.2) can be solved to yield an equation analogous to (1.15), 


a+ = A(z)P+, (2.3) 

where A is R x l'*'. Note that although wave modes on the left of the 
interface do not appear in (2.3), the projection process by which this 
equation is obtained Imposes the condition that the wave energy on the left is 
nonzero in the leftgoing components only. In other words (2.2) - (2.3) 
describes the response of the Interface to energy incident on the right. 

Now consider a finite difference model Q with p interfaces located at 
grid points j = jp**»,jp, and write Aj = j^ - j^. (To be precise, each 
^ half —integer , with one difference formula applied on j j ^ j ^ j 
and another on j^ < j < j e z. ) in this section the indices j^ are 

to be kept fixed as h,k ^ 0, and we recognize this assumption by calling ^ 
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a modal of '*fixed—Aj" type* As mentioned in the Introduction, a fixed— Aj 
problem might come up in the analysis of adaptive mesh-refinement procedures, 
or with any boundary or interface discretization that involves more than two 
distinct difference formulas. We obtain the following stability criterion; 

PROPOSITION 2. A fixed- A.j multi- interface difference model is unstable 
if and only if for some |z| > 1 it admits a nonzero solution z <{)j 
containing only lefteoing waves to the left and rightgoing waves to the right 
of all interfaces . 

Proof . The situation is illustrated in Fig. 3. For a proof, one can 

relabel the grid points so ^hat the interval from to jp becomes one 

complicated interface separating the two regions j < Ji and j > jp* Then 
the folding argument mentioned above for a single interface applies. || 


/I 





(a) interior (b) boundary 

Figure 3. Unstable multi- interface solutions z’^ 4»j at 
an interior interface and at a boundary. 


Remark . In the case of an initial boundary value problem with a boundary 
at the left, say, the region to the left of the Interfaces in Prop. 2 becomes 
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finite in extent (or possible empty, depending on labeling), so in principle 
one should not restrict the search for unstable modes to solutions that are 
leftgoing there. But in this region the difference formula is necessarily 
one-sided, which implies under the usual assumptions that it admits leftgoing 
waves only for |z| ^ 1. Therefore the change is vacuous. 

From the wave propagation point of view the following result should now be 
unsurprising. 


PROPOSITION 3. j;or the stability of a fixed- Aj multi- interface model. 
^^fff^isnt that the individual interfaces be stable . 


Remark. Stability of the Individual interfaces is presumably not 
necessary, either. 

Proof^. The proof consists of exhibiting Example 2, below, but the idea 
behind it is Indicated in Fig. 4. Imagine two interfaces at which waves can 
reflect with a reflection coefficient greater than 1. When these are placed 
together, it might happen that the reflected wave from each Interface serves 
to stimulate the reflected wave from the other. A process of reflection back 
and forth will then ensue in which at each circuit, the amplitude grows by a 
factor const > 1. Since one circuit takes only a fixed number of time steps, 
this process will cause growth at a rate Hv^ll = const’^, which is an 
explosive instability. m 


- 18 - 


(a) 


(b) 


Figure 4. The concatenation of stable interfaces may be unstable. 


Example 2 . Let = u^ on x > 0 be modeled by an "interior" formula 
for j > 2 combined with additional boundary formulas at j = 0 and j = !• 
The interior formula is an upstream difference with some added dissipation; 


v”« . v5 + - vj) + ^ - 2V" + v^.,) J . 2. (2.«) 


The formula at j = 0 is a linear combination of upstream differences: 


n n 


- v" 


n+1 n . X r 2 ^>) + — [ ^ • 


n+i n , A 
^0 = ^0 + 8 ^ 2 


(2.5) 


At j = 1 we use a leapfrog formula with some added dissipation: 

vn+l ^ ^n-1 + ^(^n _ ^n) ^ ^(^n+l _ 2^^+! + (2.6) 

It is verified in Sec. 6.3 of [24] that if ^ =J and e = 1036/83205, then 
(2.4) - (2.6) is exponentially unstable, admitting a solution = z <t>j 

with z = 129/128. The eigensolution <\> has the form (1/4, 1, 1/2, 1/4, 1/8, 
•••), and can be viewed as the superposition of leftgoing and rightgoing waves 
represented in Fig. 3b. A numerical experiment confirms that (2.4) — (2.6) is 
highly unstable and blows up like (129/128)’^ [24]. // 
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We have chosen such an unwieldy example because it is contrived to have a 


special additional property: all of the formulas (2.4) - (2.6) are 

dissipative . This is of interest because as a matter of practical experience, 
dlssipativity often ensures stability. For the case of a single interface it 
has been proved by Ciment [6] (for interfaces) and later by Goldberg and 
Tadmor [9] (for boundaries) that under reasonable hypotheses, this is always 
true. See also Sec. 6.2 of [24]. Later it was claimed by Oliger [17] that 
the same must hold with multiple Interfaces. However the example above shows 
this is not so. We formulate this conclusion as a new proposition: 

PROPOSITION 4. In a fixed-A j model with two or more interfaces, such 
as an init ial boundary value problem model with distinct boundary conditions 
^ j = 0 and j = 1, dlssipativity of each individual difference formula is 
not sufficient to ensure stability. 

It would of course be more satisfying to find an Illustration of this 
principle that was somewhat realistic. 

Example 2 also serves to Illustrate another (weaker) stability principle. 
In some circles, where the Kreiss/Osher theory is considered too complicated 
for practical work, the "von Neumann" or "Fourier method" for heuristic 
stability analysis is used Instead. This idea, proposed by Trapp and Ramshaw 
[23] (not by von Neumann) , is to check the numerical boundary formulas for 
amplification factors greater than 1 just as if they were interior formulas , 
and hope that if there are none such, the model will be stable. In general 

there is little reason to expect this procedure to work, and indeed the 
heuristic justification of it by Trapp and Ramshaw is not really valid. Yet 
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because of the algebraic simplicity of the difference formulas usually 
encountered, the idea is surprisingly reliable in practice [22]. In 
particular, for a dissipative difference model with a solvable boundary 
condition applied at a single point, it can readily be shown that the Fourier 
condition is sufficient for stability [9]. 

But Example 4 confirms that the same does not hold when there is more than 
one boundary condition: 

PROPOSITION 4' . In an initial boundary value problem model involving 
distinct boundary conditions at j = 0 and j = 1, the "von Neumann method" 
of boundary condition analysis is not sufficient to ensure stability . 

If the stability of each interface individually is not enough for a 
general stability test, what is? The unfortunate answer is that for a 
complete analysis one must investigate all possible modes z <|)j suggested by 
Props. 1 or 2 to see if they satisfy the boundary conditions. The difficulty 
with the computation is that its size grows with the total width of the 
Interface region: one must study a matrix function E(z) of dimension 
approximately Aj in the scalar case, NAj in general. The required 
investigation can be prohibitively difficult. 

However, various sufficient but not necessary conditions for stability can 
be derived that involve the interfaces individually. Consider the two- 
interface model Q = Q_lQolQ+ illustrated in Fig. 5. Here Q_, Qq, and Q+ 
are constant— coefficient difference formulas with stencil parameters » 
{Jl,r}, {Jl+,r^}, , and and j 2 (> Ji) are half-integers with 

= interface at consists simply of an abrupt change 
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from Q_ applied for j < 
^ 2 * 


to 


Qq applied for j > and similarly at 



Figure 5. Two— interface model (fixed Aj) 


We assume that each interface Q_|Qq and QolQ+ is individually stable, and 
seek a condition to ensure that no imstable solution (1)^ of Q with 
|z| > 1 (as in Fig. 4b) can exist. To ensure decoupling of Q_ and Q^, we 
assume further r_ < r + Aj and < A + A j . 

Let Qq admit R rightgoing and L leftgoing solutions, labeled as in 
Let K(z) and M(z) be the R^R and L^L nonsingular matrices 

K(z) = diag(Kj,**«,iCj^), M(z) = diag(pj,**»,Pj^), 

and let V and W be the N ^ R and N ^ L matrices with columns Vjjj and 
Wm» 

V(z) = W(z) = (Wj,“*,Wl). 

Then (1.12) can be rewritten 
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z”[v(z)K(z)J a + W(z)M(z)^ p]. 


(2.7) 


n 


By definition of V, K, W, and M, this expression satisfies Qq for all j , 
regardless of a and p. Conversely, a function = z^ <t>j satisfies 

Qq for < j < only if it has a representation (2.7) valid in 
jj^-J!,<j<j 2 +r for some a and p. The question is, for which a and 
P, if any, can a function Vj defined by (2.7) in 

extended to a solution of Q for all j that is leftgoing in j < and 

rightgoing in j > j 2 ? The answer is; For precisely those a,P satisfying 
the reflection equations 


a = Aj B, p = A 2 a (2.8) 

where Aj^ is an R^L matrix as in (2.3) relating a to P at the Q_IQq 
interface, and A 2 is an L^R matrix relating P to a at Qq|Q^_. This 

follows from the construction of (2.2). The assumption that each interface is 
stable in Isolation has permitted us to pass from the form (2.2) to (2.3), 
since it implies by Prop. 1 that Ej(z) and £ 2 ( 2 ) are nonsingular for each 
|z| > 1. 

The matrix A(z) of (1.15) was effectively defined with respect to the 
grid point j = 0, in the sense that it is at that point where a solution 

(1.12) to Q has the form Va + WP with a = A(z)P. For the present 

problem, it is more natural for Aj to be defined with respect to the grid 
point j^, and A 2 with respect to ^2’ accomplish this by 

~ j 1 J ~ j 2 •^2 

replacing A]^ in (2.8) by K Aj^ M and A 2 by M A 2 K . Equation 
(2.8) becomes 
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(2.9) 


h Ji 

K(z) a = Aj(z)M(z) ^ p, 


22 Jo 

M(z) P = A 2 (z)K(z) a. 


( 2 . 10 ) 


With this somewhat cumbersome notation it is possible to state a simple 
lemma on the existence of solutions (J)^ to Q. 


LEMMA 1. The fixed- Aj two-interface model Q described above admits a 
solution <t>j with | z | > 1 consisting of outgoing waves only in 

j < jj and j > J2 if and only if the L><L matrix 


Ej^(z) = M(z)"^j A 2 (z)K(z)^j A^ (z) 


( 2 . 11 ) 


has an eigenvalue !• 


Proof , Suppose Q has a solution Vj = ^ of the kind described. 

Let a and p be the coefficient vectors for the representation (2,7) of 

V Py definition of and A 2 > the equations 

(2*9) and (2,10) must hold. Multiplying them together gives 


M(z)^^ p = A2 (z)K(z)^^ A^(z)M(z)^^ P, 


that is , 


[m(z)^^ P] = Ej^(z)[m(z)^^ P]. 


Thus M(z) P is an eigenvector of the sort required. 
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times a 


Conversely, if eigenvalue 1, let P be M(z) 

corresponding eigenvector, and define cc by (2»9)» Then by definition of P, 
(2.10) is satisfied also, so Q has a solution of the required kind. B 

Lemma 1 now makes it possible to give sufficient conditions for stability 

based on Aj^ and A£ alone. 

PROPOSITION 5. In the fixed- A.1 two-interface problem described above, in 
which each interface individually is stable, a sufficient condition for 
stability is 

HAj^(z)II < 1, HA 2 (z)H < l for all |z| > 1 

in any norm II • II subordinate to a vector norm, if at least one of the two 
inequalities is strict for each z . 

Remark . Aj^ and A 2 are rectangular matrices, l.e., operators 

and A 2 :C^ ^ The norms in Prop. 5 are the operator norms 

subordinate to any norms on and C^, which must however be the same for 

both A 2 and A 2 « 

Proof . By the definitions of rightgoing and leftgoing we have 
Ik I < 1 < |u I for all z and m, hence IlK(z)II, llM(z)“^H < 1 in any norm. 
Together with the hypotheses and (2*11) this implies llEj^(z)H < 1 for each 
|z| > 1, which precludes the existence of the eigenvalue 1 of Lemma !• B 

Example 3 * Here we reproduce a ”P-stability” result of Beam, Warming and 
Yee [1] by considering reflection coefficients. Let u^ = u^^ on [0,1] be 
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modeled by any of the "A-stable" foraulas Q of Beam and Warming, which 
consist of the usual three-point difference operator in x coupled with an A- 
stable linear multistep formula in t. Examples are the backward Euler and 
trapezoidal (= Crank -Nicolson) formulas. Let the boundary conditions be 
= 0 at X = 1, j = Aj + 1 > 2, and q^^-order space extrapolation 
(q < Aj +1) 

(K-l)*! = 0 (2.12) 

at X = j = 0, where K denotes the shift operator Kv') = v*?. , . We claim 

J J+1 

that for any fixed Aj, Q admits no solution (])j with |z| > 1. 

Since the spatial difference in Q is just (K - K“^ , it is readily seen 

that for each |z| > 1, Q admits one rightgoing wave z^ and one 

leftgoing wave , with Re< < 0 < ReH, |<| < 1 < ||i.|, and = -l/<. 

The first inequality is derived as follows in Theorem 2.4.1 of [24]. If Q 

is A-stable, then Re(K-l/ic) < 0 implies |z| < 1. Contraposltively , |z| > 
1 implies Re(K - 1 /k) > 0. Since |k| < 1, this means |z| >1 implies 
Rex < 0. 

Now we compute reflection coefficients. At j = A j+ V 2 one has 

A 2 = - / K/p = - ix, (2.13) 

and at j = V 2 » 

Aj = - / ( ^ Z — • (2.14) 

1 “4" 1C 

By the above inequalities one has | _ ^ | < 1 for |z| > 1, and therefore 
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For q 


IA 2 I < 1, lAj < \K\^-^ . 

= 1 both reflection coefficients have magnitude < 1, and by the 

argument of Prop. 5 we are done. If q > 1, the assumption Aj + 1 > q 

Ai 

implies that the term K in (2.11) cancels any amplification due to the 
factor l<|^ ^ above, so stability follows from Lemma 1. Alternatively, to 
stick with the one-boundary-at-a-time approach of Prop. 5, one can renumber 
the vertices so that the lefthand boundary lies at j = q- V 2 instead of j = 
V 2 > and then lAj^l will be <1 regardless of q. // 

Remark. A similar argument can be applied to the LF model (1.8) 
together with a space-time extrapolation condition such as ~ ^1* 


3. TWO INTERFACES SEPARATED BY A FIXED DISTANCE Ax 

In this section we continue to investigate the configuration illustrated 
in Fig. 5, except that Ax rather than Aj will be held constant. Consider 
a two— interface model Q = Q_ IQq I Q_|_ in which the interfaces lie at 
positions Xj^ = j ^ h and X 2 = J 2 ^ad set Ax = X 2 - Either or both 

of the interfaces may in fact be a boundary; if both of them are, then Q is 
a model for a differential equation on a strip such as [0,1]. We ask: as the 
mesh is refined, i.e., as h,k 0 with and X 2 fixed, will the 

behavior of Q be stable or unstable? 

It is now that the distinction between stability and P-stabllity becomes 

important. Following Beam, Warming, and Yee [1], define: 
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Definition « The fixed-Ax two-interface model Q described above is 

P-stable if it is GKS-stable and in addition, for each fixed h > 0 it admits 
no solutions Vj = with |z| >1 containing only leftgoing waves to 

the left and rightgoing waves to the right of both interfaces. 


(”P” stands for "practical".) Actually, P-stability is not a stability 
concept of the usual sort, since it is defined in terms of what eigensolutions 
Q admits rather than what growth estimate it satisfies. But obviously this 
condition is vital if the time-dependent finite-difference model is to be used 
to approximate steady-state solutions, a procedure that is common in practice. 
In their tests Beam, et al. found P-stability of a linearized model problem to 
be a good indicator of success in practical nonlinear steady-state flow 
calculations [31] . 

We begin with the following result due to Kreiss: 

PROPOSITION 6. The flxed- Ax two-interface model described above is GKS- 
stable if and only if both interfaces Q |Qq and QqIQ^ are individually 
GKS-stable. 

Proof . See Section 11 of [10] and also Section 2 of [12]. The result 
refers specifically to GKS-stability , and is not necessarily valid for other 
definitions such as Jl 2 ”Stability . The basis of the argument is the 

invariance of GKS-stability with respect to perturbations of size 0(k); the 
effect of each boundary on the other can be shown to be of this order 


The conclusion of Prop* 6 corresponds to vhat is often observed in 
practice: if each of two interfaces is GKS-stable, the computational results 

are usually satisfactory, while if one of them is not, they are usually wrong 
and sometimes explosively so* But this section can be viewed as an 
investigation of how Prop* 6 fails to tell the whole story* Our remaining 
results can be summarized as follows* Proposition 7 shows that repeated 
reflections between GKS-stable interfaces can cause P-unstable growth at a 
rate const^, even though GKS-stability is maintained (cf* [1] and Section 7 of 

[10] )* Proposition 8 shows that reflection between weakly GKS-unstable 

, , , const t , _ 

interfaces can cause catastrophic growth at the rate (ct* 

Section 17 of [14])* Proposition 9 shows that the latter problem will not 

occur when the unstable interfaces have finite reflection coefficients* 

Proposition 10 shows that dissipation prevents both growth phenomena (cf* 

[11] )* Finally Prop* 11, like Prop. 5, shows that all growth can be ruled out 
if the reflection matrices at the two interfaces are known to satisfy 

IIA^II < 1* 

PROPOSITION 7* GKS-stability does not imply P-stability* Specifically, 
let each Interface in the flxed- Ax two-interface model described above be 
GKS-stable* If the reflection matrix at one or both interfaces has norm 

greater than 1, then repeated reflections between the interfaces may lead to 
solution growth at the rate 

llv“ll > (const)*^ llv°II. (3.1) 

Proof . In the following discussion we first explain the growth rate 

const^ by two different heuristic arguments, which will be used again later in 
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this section. The purpose of these arguments is to show that, although growth 
at the rate (3.1) need not occur for every model satisfying the hypotheses, it 
is nevertheless the typical growth rate to be expected in such problems. The 
proof of the proposition as stated then consists of exhibiting Example 4. H 


Argument by repeated reflections . The principle of Prop. 7 is the same 
as that of Fig. 4, except tht Ax rather than Aj is held constant. Suppose 
that for some |z| = 1, a (nondissipating) wave of frequency z exists which 
can travel leftwards with C < 0, reflect at the Q_|Qq interface into a 
rightgoing wave with C > 0, and then reflect at the QqIQ+ interface into 
the original leftgoing wave mode again. If the product of the two reflection 
coefficients in this circuit is greater than 1, then amplification has taken 
place, and it will be repeated in further reflections. The time taken to 

complete each circuit is roughly constant, independent of h and k as 
h,k 0. Therefore one must expect growth at the rate const^« 


Argument by perturbed reflection coefficients . If Q permits geometric 
growth in t, we can expect the existence of an eigensolution Vj = z^ (t)^ 
with |zq| > 1; the rate of growth will depend on how large |zq| can be. 
For simplicity suppose that Qq admits just one rightgoing mode and 

one leftgoing mode z^ \x^ for each |z| > 1, and as in the above argument, 
suppose that for some |zq| = 1 one has \<\ = ||i|=l, and 

lAf A^l > !• Then the diagonal matrices K and M of Section 2 reduce to 
K and p, and the matrix of (2.11) is a scalar with modulus jA^ A 2 I. 

Obviously this scalar is not equal to 1, so by Lemma 1, Q does not ,have a 
solution Zq 4>j • But suppose it happens that Ej^ = 1 + const, where, here and 
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from now on, const denotes a quantity of order of magnitude 1 that varies 
from one occurrence to the next and is positive except possibly for an 
imaginary part of size 0(e), when this makes sense in context. To find a 
solution satisfying (2.7), consider z = Zq( 1 + e), 0 < e « 1. This 

perturbation changes k, \i, and A 2 by 0(e). In particular k and \i 
become 

K ->• k( 1 - conste), p p(l + conste) . 


(In the limit e = 0 the constants here are lAICj^l and 1A|C^|.) By 
(2.11) therefore becomes 


E = (1 + const) (1 - conste) 

Ij 


Aj 


For El to have value 1, the two factors have to balance, which means 
e = 0(1/Aj). Therefore one can expect that any eigensolution <l>j to Q 

will grow at the rate 


|z|' 


U + 


Aj 




as asserted in (3.1). 


Example 4 . Let u^ = u^j on [0,1] be approximated by the LF formula 
(1.8) together with the (admittedly contrived) boundary conditions 


V 


n+1 

0 



+ v! 


n-2 




Aj+1 


= 0 . 


(3.2) 
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The reflection coefficient functions are easily seen to be 


A 


1 


- '^< 7 ^ 






(3.3) 


and since the denominators are never zero, both interfaces are GKS-stable. 
However, lAjI can be larger than 1. For simplicity consider the semidiscrete 
limit X = 0, z = 1. By (1.17), for any 6 = Ch e [0,ii/2), LF then has a 
solution 


z = 1, 



C = -COS0, C = COS0* 


For any 0 with |Aj^(0)| > 1, one can expect Q to admit an eigensolutlon 
that grows approximately at the rate | A^ (0) | since each circuit of a 
trapped wave packet will take time 2/cos 0. The maximum of these rates for 
the given formulas turns out to be at 0 « .75, where one gets |Aj| « 2.38, 

C " .725, and growth (1.37)^. Numerical experiments confirm that solutions 
grow roughly at this rate, independent of h and k. 

To establish Prop. 7 rigorously, one must prove that Q admits the kind 
of growing eigensolutlon we have described. This can be done by using pertur- 
bation arguments based on the above heuristic reasoning to show that (2.11) in 
Lemma 1 has a solution with |z| « (1.37) • Since the conclusion is so 
obvious, we will not give details. // 

The possibility of P-unstable growth as in Example 4 was recognized from 
the start by Gustafsson, Kreiss, and Sundstr'om, and in fact Section 7 of [10] 
is devoted to determining \dien it will occur in a certain 2><2 problem.. In our 
particular example, the model remains P-unstable no matter how small h and 
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k become. Beam, et al. give the Impression in various papers that this 
cannot happen, but that is true only when one is dealing with dissipative 
formulas; see Prop. 10 below. The reason that dissipation did not ensue 
P-stability for the values of h and k they were dealing with was that, 
because of their interest in steady-state results, they were using very large 
values of and their formulas happened to be nondissipative in the limit 
X Thus their computations made use of difference formulas that were 

dissipative but only weakly so. 

Now let reflection coefficients be present that are not merely greater 
than 1, but infinite. The potential growth rate becomes much more severe. 

PROPOSITION 8- Let one or both interfaces in the fixed- Ax two-interface 
problem be algebraically GKS-uns table, with an infinite reflection 
coefficient. Then repeated reflections betwen the interfaces may sometimes 
lead to solution growth at the exponential rate 

Hv"ll > ^ llv°ll. (3.4) 


Remark . For a single GKS-unstable interface with an infinite reflection 
coefficient, it is shown in [27] that the unstable growth is in general no 
worse than ilv^II « const n liv^ll. This is what is meant above by 

"algebraically" GKS-unstable. 

Proof . Again we will motivate (3.4) by two arguments. Then we prove the 
proposition by exhibiting Example 5. |H 

Argument by repeated reflections * Suppose Aj^ is infinite at z = zq, 
and behaves near there like 
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(3.5) 


Since there are only Aj points between the interfaces for each fixed h, 
Fourier analysis implies that no wave on the Qq grid can have a spectrum 
narrower than 0(1/Aj), Therefore it is plausible that in applying (3.5) to 
the finite grid, the largest amplification possible will be that obtained with 
an effective value with ~ I = const/Aj, i.e., HA^II *> const A j . 
Since as before each circuit takes roughly a fixed amount of time, independent 
of h and k as h,k 0, this leads immediately to (3.4). 


Argument by perturbed reflection coefficients . As before, suppose that 

Qq admits one leftgoing mode z“ and one rightgoing mode z“ for each 

|z| > 1, and that for some |zq| = 1 one has |K| = |p,|=i^C <0<C, 

|i K 

IAiCzq)! = “, and |A2 (zq) ( > 0. Suppose furthermore that Aj^ behaves like 
(3.5) for z « Zq. Then under the perturbation z = Zq( 1 + e) one has 

A • 

K ->■ K(1 - conste), p ^ p(l + conste), E ~ ^°’^ste) ^ ^ 

L e 

For = 1 we must have (1 - conste)^J = e, which Implies 


^ ^ const ^ / A • \ 

e « rr— log(Aj) 


Aj 


Hence growth should be expected at the rate 


llv^ll 


= Z 


|n ^ ^ constlogAj ^ tA.l ^ ^tconstlogAj ^ ^^.^constt 


Aj 
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Example 5 ([14], Sec. 17). Let u^. = on [0,1] be approximated by 

LF (1.8) with boundary conditions ~ (1*16) and ~ 

have seen in (1.18) that this model has an Infinite reflection coefficient 
at z = -1 ; in fact one has as in Example 3 


Aj(z) = -/k/p I 3 k ’ 


= -/k/P . 


With these formulas (2.11) becomes 


Ej^(z) = (K/p)^j+^ j- 




- 1C 


and since p. = -1 /k for LF, this can be rewritten 


E^iz) = (-K^)W 


Assume Aj+1 is even, and write ic = 1 - 6. The condition Ej^(z) = 1 

becomes 


(1 - 6 ) 


2Aj+2 _ 6-6^ 


2 - 6 * 


It is obvious that this equation has a positive real solution near 6 = 0, 
which is asymptotic to 6 = log Aj/2Aj as Aj “. The corresponding value 
of z is asymptotic to 




Therefore Q has an eigensolution which grows at the rate 
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(1 . (l ^ 10|^)tAjA . ^tlogij/2 ^ 

This matches the result stated as (17.10) in [14], and numerical experiments 
confirm that physical solutions are rapidly obliterated by growth at the 
predicted rate* // 


The possibility of catastrophic two-boundary interactions as in Prop* 8 
has long been recognized by Kreiss and his colleagues, and it has been given 
sometimes as a justification of the apparent strictness of the GKS stability 
definition* We now show that this justification is only partial, for not all 
GKS-unstable boundaries have infinite reflection coefficients, yet an infinite 
reflection coefficient is required for the catastrophic two-boundary 
interaction to occur: 


PROPOSITION 9* Let one or both interfaces in the fIxed- Ax two -in ter face 
problem be algebraically GKS-unstable, but with finite reflection coefficients 
only* 


(const)^. 


Then Q admits no eigensolutions Vj = z <t)j that grow faster than 


n 


Proof * Consider an eigensolution z <J)j of Q, and let M, K, A 2 , be 

the matrices of Lemma 1 for the given value z* By Lemma 1, the matrix 




A2 


(3.6) 


has an eigenvalue 1, which implies >1. On the other hand, the 

finite reflection coefficients assumption implies 
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for some T < “ 


These bounds together yield 


HM ^11^^ hkh^^ > ^ 

oo oo T 

or in particular, since |< | < 1 < 1 m| for each of the entries in M and K, 

J_ 

l>=l. |u‘‘l > (3.8) 

for some k and p. 

Now the critical observation is that for any Cauchy stable formula Q, |z| 
-1 is bounded by a multiple of 1 - |k | when the latter is small. For a 
proof, see Lemma 9,1 of [10]; the constant factor is essentially X times the 
maximum group velocity admitted by Q, Therefore the last inequality implies 

^ ^const/Aj 

for large enough Aj. But this leads to 

n const 

I z 1^ < (T) = const^, 

which proves the theorem, 

Now we come to the question of dissipation. In the fixed-Aj situation, 
the use of dissipative formulas gave no guarantee of stability, because the 



attenuation introduced by dissipation might always be overcome by 
amplification due to reflection at the boundaries. But in the fixed-Ax 
problem, the attenuation of any nonphysical wave mode will Increase as the 
mesh is refined. For fine enough meshes, this must overcome any finite 
amplification factors. 

PROPOSITION 10. Let Q be a GKS-stable, totally dissipative, consistent 
model of a fixed-A x two-interface problem. Let the solutions to the problem 
being modeled satisfy a bound 

llu(t) II < const e ^ Hu (0)11. 

Then for all sufficiently small h and k , Q is P-stable . 

Remark. By "totally dissipative," we mean that the interior model Q 
dissipates oscillations with respect to t as well as x. For two-level 
formulas this is the same as the usual definition of dissipativity . For 
multilevel formulas, there is the additional requirement that the scheme admit 
no solutions v^ = 4), (J) = const, with |z| = 1 but z ^ 1 [24]. 


Remark. This is an elaboration of the theorem stated by Gustafsson in 
[11]. See also [7]. 


Proof . We must show that Q admits no eigensolution 4)^ with 

|z| > 1, for large enough Aj. Suppose to the contrary that for a sequence of 
values Aj -► <=, Q has a solution z“ 4>j with |z| > 1. Since Q is GKS- 
stable, it has finite reflection coefficients, so (3.6) - (3.8) of the last 
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proof are again valid* Equation (3*8) Implies |k| ^ 1 and |p»| 1 as 

Aj -> ®* By dissipativity , this implies k 1 and |i -> 1. By total 
dissipativity, this in turn implies z ■> 1 also. In other words, as 
Aj -► any eigensolutions of Q have to approach physically meaningful 
solutions with low frequency and wave number. 

Let the matrix A in the differential equation (1.1) have L positive 
and R = N - L negative eigenvalues. (There can be no zero eigenvalues, or 
the decay assumption would fail.) Consider the behavior of the L^L matrix 
of (3.6) as Aj ® and z ->• 1. By consistency, L values and R 

values K approach 1, while the remainder are bounded away from 1 in 

modulus. The powers and in (3.6) therefore approach diagonal 

/v» ^ 

matrices containing L and R ones and L — L and R — R zeros, 
respectively. On the other hand, the reflection matrices Aj and A£ also 
approach limiting values, and by consistency, the restrictions of these to the 
L + R rows and columns corresponding to modes k = 1 or P = 1 must equal 
the R^L and LXR reflection matrices Aj^ and for the differential 

equation being modeled. Combining these observations, we obtain 


^ as Aj ■> ”, 


where E^ is an L^L matrix consisting of the L^L matrix A^ A 2 padded with 
additional rows and columns of zeros. 

Now by the decay assumption, the eigenvalues of must be bounded below 

1 in modulus; otherwise some eigenvector would describe a solution to the 
differential equation consisting of a wave that reflected back and forth 


without decaying. On the other hand by Lemma 1, we have assumed Ej^ 
eigenvalue 1. This is a contradiction. 


has an 
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Our final result Is the same as Prop. 5, but restated for the fixed- Ax 
problem. 


PROPOSITION 11. In the fixed- Ax two-interface problem, in which each 

interface individ ually is GKS-stable, a sufficient condition for P-stabilitv 

is 

IlAj (z) II < 1, IIA 2 (z)II < 1 for all |z| > 1 

is any no rm subordinate to a vector norm, if at least one of the two 

inequalities is strict for each z. 

Proof . Same as for Prop. 5. HI 
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